#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <complex.h>
#include <string.h>
#include <unistd.h>

#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_vector_complex.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_fft_complex.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_sort_vector.h>

#define N 800


// gsl_vector_view gsl_matrix_subrow (gsl_matrix * M, size_t I, size_t OFFSET, size_t N)
// size_t gsl_vector_max_index (const gsl_vector * V)

//void histograma(gsl_vector *eval, gsl_matrix *hist, int n);

void obterDiscos(gsl_matrix* A, gsl_matrix* info);
void histograma(gsl_matrix* info, gsl_vector* hist, int m);
void obtersol(gsl_vector* histGeral, gsl_vector* hist, int m);
void gerarNovaMatriz(gsl_matrix* A);
int numAutovals(gsl_vector* histGeral, int m);
int histEigen(gsl_vector* hist, gsl_vector* autoval);

void imprimirHistograma(FILE *fp, gsl_vector* hist, int m);
void imprimirMatriz(gsl_matrix* A);
void printVetor(gsl_vector* autoval, int num);
void imprimirRaios(gsl_matrix* info);

int main(void)
{
	int i, j, k, l, lk;
	int i1, j1;
	gsl_matrix *A, *info, *H;
	gsl_vector *hist, *histGeral;
	gsl_vector* autoval;
	double tmp, theta;
	FILE *fp;
	char str[100];
	int m, num = 0;
	gsl_eigen_symm_workspace * w;
	double delta = 0.000001;
	int rd;
	int vezes = 0;
	
	m = 20000;

 	srand(time(0));
	
	A = gsl_matrix_alloc(N, N);
	gsl_matrix_set_zero(A);
	H = gsl_matrix_alloc(N, N);
	hist = gsl_vector_alloc(m);
	histGeral = gsl_vector_alloc(m);
	info = gsl_matrix_alloc(N,2);
	
	for ( i = 0 ; i < m ; i++ )
	{
		gsl_vector_set(histGeral, i, N);
	}
	

	gsl_matrix_set(A, 0, 0, 1);
	gsl_matrix_set(A, 0, 1, 2);
	gsl_matrix_set(A, 1, 0, 2);
	for ( i = 1 ; i < (N-1) ; i++ )
	{
		gsl_matrix_set(A, i, i, (2*i+1)/600.0);
		gsl_matrix_set(A, i, i+1, (i+4)/600.0);
	}
	gsl_matrix_set(A, N-1, N-1, (2*N-1)/600.0);
	for ( i = 1 ; i < (N-1) ; i++ )
	{
		gsl_matrix_set(A, i+1, i, gsl_matrix_get(A, i, i+1));
	}

	w = gsl_eigen_symm_alloc(N);
	autoval = gsl_vector_alloc(N);
	gsl_matrix_memcpy(H, A);
	gsl_eigen_symm(H, autoval, w);
	gsl_sort_vector (autoval);
//	printVetor(autoval, N);

	do
	{
		vezes++;
		gerarNovaMatriz(A);
		obterDiscos(A, info);
		histograma(info, hist, m);
		if ( num == 0 )
		{
			for ( l = 0 ; l < m ; l++ )
			{
				gsl_vector_set(histGeral, l, gsl_vector_get(hist, l));
			}
		}
		obtersol(histGeral, hist, m);
		num = numAutovals(histGeral, m);
		printf("num = %d\n", num);
		if ( vezes == 20 ) break;
	} while (num != N );
	
	sprintf(str, "resultado/Histograma-teste", num);
	fp = fopen(str, "w+");
	imprimirHistograma(fp, histGeral, m);
	fclose(fp);
	
	printf("RAIOS\n");
	imprimirRaios(info);
	
//	printf("AUTOVALS\n");
//	num = histEigen(histGeral, autoval);
//	printVetor(autoval, num);
	
	return 0;
}

int histEigen(gsl_vector* hist, gsl_vector* autoval)
{
	int i,j,k,m;
	
	m = 20000;
	
	for ( i = 0, j = 0 ; i < m ; i++ )
	{
		if ( gsl_vector_get(hist,i) != 0 )
		{
			for ( k = 0 ; k < gsl_vector_get(hist,i) ; k++ )
			{
				gsl_vector_set(autoval, j, (i-m/2)/100.0);
				j++;
			}
		}
	}
	return j;
}

int numAutovals(gsl_vector* histGeral, int m)
{
	int num = 0;
	int i;
	
	for ( i = 0 ; i < m ; i++ )
	{
		num = num + gsl_vector_get(histGeral, i);
	}
	return num;
}


void obtersol(gsl_vector* histGeral, gsl_vector* hist, int m)
{
	int i;
	for ( i = 0 ; i < m ; i++ )
	{
		if ( gsl_vector_get(hist, i) < gsl_vector_get(histGeral, i) )
		{
			gsl_vector_set(histGeral, i, gsl_vector_get(hist, i));
		}
	}
}

void gerarNovaMatriz(gsl_matrix* A)
{
	int i,j,k;
	double s, c, eps, delta;
	double Aii, Aij, Ajj, Aki, Akj, ss, cc, sc, theta;
	int rd, itmp;
	double tmp;
	static int vezes = 0;
	
	eps = 0.0001;
	
	delta = 0.1/pow(1.5,vezes);
	vezes++;
	if ( vezes > 20 )
		vezes = 20;

	
	if ( vezes == 1 )
		rd = rand()%(N-1);
	else
		rd = 0;
	
	for ( i = rd ; i < N ; i++ )
	{
		/* Qual j ? */
	//	tmp = -1;
		for ( j = i+1 ; j < N ; j++ )
		{
			
// 			if ( fabs(gsl_matrix_get(A, i, j)) > tmp )
// 			{
// 				tmp = fabs(gsl_matrix_get(A, i, j));
// 				itmp = j;
// 			}
// 		}
// 		j = itmp;
			Aij = gsl_matrix_get(A, i, j);
			if ( fabs(Aij) < delta )
				continue;
// 	
			Aii = gsl_matrix_get(A, i, i);
			Ajj = gsl_matrix_get(A, j, j);
			
			if ( (fabs(Aii-Ajj)) >= eps )
			{
				theta = 0.5*atan(2.0*Aij/(Aii-Ajj));
			}
			if ( (fabs(Aii-Ajj)) < eps )
			{
				theta = (M_PI/4.0)*(Aij/(fabs(Aij)));
			}
			if ( (fabs(theta)) > M_PI/4.0 )
			{
				theta = theta - (M_PI/2.0)*(theta/(fabs(theta)));
			}
			//printf("theta = %f\n", theta);
			c = cos(theta);
			s = sin(theta);
			ss = s*s;
			cc = c*c;
			sc = s*c;
			gsl_matrix_set(A, i, i, cc*Aii + 2*sc*Aij + ss*Ajj);
			gsl_matrix_set(A, j, j, ss*Aii - 2*sc*Aij + cc*Ajj);
			for ( k = 0 ; k < N ; k++ )
			{
				if ( (k == i) || ( k == j) ) continue;
				Aki = gsl_matrix_get(A, k, i);
				Akj = gsl_matrix_get(A, k, j);
				gsl_matrix_set(A, k, i, Aki*c + Akj*s);
				gsl_matrix_set(A, i, k, gsl_matrix_get(A, k, i));
				gsl_matrix_set(A, k, j, -Aki*s + Akj*c);
				gsl_matrix_set(A, j, k, gsl_matrix_get(A, k, j));
			}
			gsl_matrix_set(A, i, j, 0);
			gsl_matrix_set(A, j, i, 0);
		}
	}
	
	
//	i++;
//	if ( i == N ) i = 0;
//	printf("i = %d\n", i);
	//	imprimirMatriz(A);
	//	sleep(2);
	return;
}


void histograma(gsl_matrix* info, gsl_vector* hist, int m)
{
	double max, min, delta;
	double infomax, infomin;
	double tmp;
	double interval, points;
	double val[m];
	int i, j, npoints, offset;
	double eps = 1e-6;
	int ind_infomax, ind_infomin;
	
/*	min = -100.0;
	max = 100.0;
	delta = 0.01;
	npoints = (max-min)/delta;
	offset = fabs(min)/delta;*/
	offset = m/2;

	gsl_vector_set_zero(hist);
	
	// -100 corresponde ao índice 0, logo +100 não faz parte do vetor
	
	for ( i = 0 ; i < N ; i++ )
	{
		infomax = gsl_matrix_get(info, i, 0) + gsl_matrix_get(info, i, 1);
		infomin = gsl_matrix_get(info, i, 0) - gsl_matrix_get(info, i, 1);
	
		ind_infomax = (int)(infomax*100) + offset;
		ind_infomin = (int)(infomin*100) + offset;
	
		for ( j = ind_infomin ; j <= ind_infomax ; j++ )
		{
			gsl_vector_set(hist, j, gsl_vector_get(hist, j) + 1);
		}
	}
}

void obterDiscos(gsl_matrix* A, gsl_matrix* info)
{
	int i,j;
	double cont;
//	printf("info\n");
	for (i = 0 ; i < N ; i++ )
	{
		gsl_matrix_set(info, i, 0, gsl_matrix_get(A, i, i));
		cont = 0;
		for ( j = 0 ; j < N ; j++ )
		{
			if (j == i) continue;
			cont = cont + fabs(gsl_matrix_get(A, i, j));
		}
		gsl_matrix_set(info, i, 1, cont);
//		printf("centro = %14.4f e raio = %14.4f\n", gsl_matrix_get(info, i, 0), gsl_matrix_get(info, i, 1));
//		sleep(2);
	}
}

void imprimirHistograma(FILE* fp, gsl_vector* hist, int m)
{
	int i;
	for ( i = 0 ; i < m ; i++ )
	{
		fprintf(fp, "%10.4f%7d\n", (i - 10000)*0.01, (int)gsl_vector_get(hist, i));
	}
}

void imprimirMatriz(gsl_matrix* A)
{
	int i,j;
	
	for ( i = 0 ; i <  N ; i++ )
	{
		for ( j = 0 ; j <  N ; j++ )
		{
			printf("%14.5lf", gsl_matrix_get(A, i, j));
		}
		printf("\n");
	}
}

void printVetor(gsl_vector* autoval, int num)
{
	int i;
	
	for ( i = 0 ; i < (num-1) ; i++ )
	{
		printf("%-4.6lf\n", gsl_vector_get(autoval,i));
	}
	printf("%-4.6lf\n", gsl_vector_get(autoval,i));
}

void imprimirRaios(gsl_matrix* info)
{
	int i;
	for (i = 0 ; i < N ; i++ )
	{
		printf("%-4.6lf\n", gsl_matrix_get(info, i, 1));
	}
}

/* gcc -lm -lgsl -O3 -lgslcblas jacobi.c -o jacobi.exe */